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Abstract 

The structure and transport properties of Si02~Al203 melts containing 13mol% and 47mol% 
AI2O3 are investigated by means of large scale molecular dynamics computer simulations. The 
interactions between the atoms are modelled by a pair potential which is a modified version of the 
one proposed by Kramer et al. [J. Am. Chem. Soc. 64, 6435 (1991)]. Fully equilibrated melts in the 
temperature range 6000 K> T > 2000 K are considered as well as glass configurations, that were 
obtained by a rapid quench from the lowest melt temperatures. Each system is simulated at two 
different densities in order to study the effect of pressure on structural and dynamic properties. We 
find that the Al atoms are, like the Si atoms, mainly four-fold coordinated by oxygen. However, 
the packing of the AIO4 tetrahedra is very different from that of the Si04 tetrahedra, which is 
reflected by the presence of triclusters (O atoms surrounded by three cations) and edge-sharing 
AIO4 tetrahedra. On larger length scales, a micro-segregation occurs, resulting in an Al-rich 
network percolating through the Si-0 network. This is reflected in a prepeak of concentration- 
concentration structure factors around 0.5 (both in the system with 47mol% and 13mol% 
AI2O3!). We also address the interplay between structure and mass transport. To this end, the 
behavior of the selfdiffusion constants for the different compositions and densities is studied. 
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I. INTRODUCTION 



One of the most abundant oxides in natural silicates and technological silicate glasses 
is AI2O3. Although alumino sihcates are therefore of central interest in geosciences and 
materials science, the structure even of the binary system Si02-Al203 is far from being 
well-understood. However, the knowledge of the chemical ordering in the latter system 
provides also the basis for a better understanding of the structure of more geologically 
relevant alkali or alkaline earth alumino-silicates. 

The chemical ordering of aluminium is very complicated when built into the tetrahedral 
Si-0 network. In order to achieve local charge neutrahty, Al^+ ions need a different envi- 
ronment of ions than Si^+ ions. Thus, unhke Si02, Al^+ ions do not form a network of 
AIO4 tetrahedra that are connected with each other via the O atoms at the corners. Instead, 
different experimental techniques such as nuclear magnetic resonance (NMR) [Lee and Steb- 
bins, 2000; Schmiicker et ai, 1999; Sen and Youngman, 2004; Stebbins et al., 2005; Stebbins 
and Xu, 1997; Xue and Kanzaki, 1999] as well as IR and Raman spectroscopy and X-ray 
scattering [Morikawa et al, 1982; Okuno et al, 2005] found evidence for structural units 
such as three-fold coordinated oxygen atoms and five- and six-fold coordinated aluminium 
atoms that are not found in pure silica, unless one considers amorphous silica at very high 
temperatures, say above 4000 K (see Horbach and Kob, 1999). 

At low AI2O3 concentrations, the Al atoms are mainly four-fold coordinated by oxygens, 
but, as proposed in NMR studies, these AIO4 units are accompanied by so-called triclusters, 
i.e. structural units where an oxygen atom is surrounded by three cations (where at least one 
of them is an Al atom) . Recently, molecular orbital calculations confirmed the possibility of 
such triclusters (Kubicki and Toplis, 2002). Moreover, a combination of molecular dynamics 
computer simulations and Hartree-Fock calculations (Tossell and Cohen, 2001; Winkler et 
al, 2004; Tossell and Horbach, 2005) found evidence that O tricluster atoms participate 
typically in two-fold rings (i.e. edge-sharing geometries) of composition AI2O2 or AlSi02. 
Certainly, further experimental studies are necessary to clarify to what extent triclusters 
and edge-sharing tetrahedra exist in real aluminium silicates. 

As shown by recent NMR studies, five- and six-fold coordinated Al atoms are also im- 
portant structural units in the system Si02-Al203. The fraction of these highly coordinated 
Al atoms tends to increase with increasing AI2O3 concentration (Sen and Youngman, 2004; 
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Stebbins et al, 2005). All these experimental findings show that, in alumino silicates, the 
local chemical ordering of Al atoms is very different from that of Si atoms. 

In a recent MD simulation study of the system (Al203)2(Si02) [abbreviated in the fol- 
lowing as AS2], we have shown that the different chemical ordering of Al and Si on local 
length scales also leads to structural ordering on intermediate length scales of the order 
of Inm (Winkler et al., 2004). This intermediate range order (IRO) can be described by 
a microphase separation where the Al-rich network structure percolates through the Si-0 
network. The IRO gives rise to a prepeak in partial static structure factors at a wave number 
q = 0.5 A~^. It can be seen as a precursor of the metastable liquid-liquid phase separation 
below ~ 1900 K that is found experimentally (MacDowell and Beall, 1969) between about 
10mol% and 50mol% (AS2, with 33mol% AI2O3, lies approximately in the center of the 
demixing region). 

The formation of IRO seems to be quite a general feature of multicomponent silicate 
melts. Similar IRO as the one found in AS2 has also been seen in neutron scattering 
experiments of a calcium sihcate (Gaskell et al, 1991), sodium sihcates (Meyer et al, 2002; 
Meyer et al, 2004), and alkali aluminosilicates (Cormier et al, 2001; Kargl and Meyer, 
2005). For the example of sodium silicates, the IRO is reflected in a prepeak around 0.9 
in static structure factors. It has its origin in the formation of sodium-rich channels in the 
static structure. These channels serve as preferential pathways in an immobile Si-0 matrix 
and thus provide an explanation for the high mobility of sodium in ion-conducting sodium 
sihcate melts (Horbach et al, 2002; Meyer et al, 2004). The presence of diffusion channels 
in sodium silicates is a nice example for the interplay between structure and mass transport 
in glassforming melts. 

The structure-transport relation is also a central issue of the present study. We have 
extended our previous simulations of AS2 to aluminium silicates with 13mol% AI2O3 and 
47mol% AI2O3. These systems are simulated at different densities to see how pressure af- 
fects the structural and dynamic properties of aluminium silicates. Furthermore, we aim to 
understand how the interplay between local structural features (triclusters, two-fold rings, 
five- and six-fold coordinated Al ions) and IRO changes with composition. This is in turn 
will allow us to elucidate the interplay between structural features and transport processes, 
when we investigate the temperature dependence of the selfdiffusion constants for the dif- 
ferent systems under consideration. 
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II. MODEL AND DETAILS OF THE SIMULATIONS 



As in our recent study of AS2, we use a potential proposed by Kramer et al. (1991) to 
model the interactions between the atoms. It is based on the so-called BKS potential (van 
Beest et al, 1990) for pure silica that has been extensively apphed in recent studies (see, 
e.g., Horbach and Kob, 1999, and references therein). Apart from mixtures of Si02 with 
AI2O3, the Kramer potential allows also the consideration of sodium silicates. As shown 
recently, it provides quite a realistic description of sodium di-, tri- and tetra-silicate (Ispas 
et al., 2002; Meyer et al., 2004) as well as of AS2 (as far as comparison with experimental 
data is possible, see also below) [Winkler et a/., 2004; Tossell and Horbach, 2005]. 

The potential has the following functional form: 

(t)ap{r) = ^"^^^ + Aaf3 exp {-Bapr) - ^ (1) 

with a, 13 = Si, Al, O. Here r is the distance between an ion of type a and an ion of type /3. 
The values of the parameters {Aa/3, B^/s, Ca/3} that were calculated by ab initio methods are 
Asio = 18003.7572 eV, Aaio = 8566.5434 eV, Aqo = 1388.7730 eV, Bsio = 4.87318 A'^, 
5aio = 4.66222 A-\ Bqo = 2.76 A-\ Csio = 133.5381 eVA^, Caio = 73.0913 eVA^, and 
Coo = 175.0 eVA6 (for the Si-Si, Si-Al and Al-Al interactions the latter parameters are all 
set to zero) [Kramer et al. (1991)]. 

The Buckingham part of the potential, 

C 

(f)Bir) = A^p exp {-Bapr) - , (2) 

has been truncated and shifted to zero at rf = 5.5 A. In order to make the truncated 
potential, 0g"'^'^(r), different iable at rf it has been multiplied by a smoothing function of 
exponential form: 



i trunc / 



r < 

(3) 

r>rf 



with cIb = 0.05 A^. The parameter rfe is chosen such that the exponential in Eq. Q does 
not affect the system's properties. But it smoothens out the cusp at which would lead 
to a discontinuous force at rf and thus to a drift in the total energy in microcanonical MD 
runs [Allen and Tildesley, 1987]. By the systematic use of smoothing functions as in Eq. Q 
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(see also below), we have not encountered an energy drift in any of the microcanonical runs 
that we have performed in this work. 

In the long-ranged Coulomb-part the charges q^e (e: charge of an electron) are not the 
bare ionic charges of ions of type a but are considered to be effective charges. The charges for 
silicon and oxygen are set to gsi = 2.4 and go = —1-2, respectively. With the original charge 
for aluminium, gAi = 1-9, the unphysical situation of a non-zero net charge results. Only 
with the additional component phosphorus the system recovers charge neutrality according 
to the parameter sets in Kramer et al. (1991). We have therefore modified the Kramer 
potential by using a distance-dependent Al charge as follows: 

= I [1 + In {Cm^^ + 1)] exp {-^) r < r,, 

with (jAi = 1-8. The concept of a "distance-dependent Al charge" should not be taken too 
literally. We just introduce an additional short-ranged potential to the original one, with a 
cut-off at r = tai. 

The exponential function in Eq. serves again as a smoothing factor. The parameters in 
Eq. (j3)) are adjusted such that the potential is very close to the original one at short distances. 
At distances r > tai the Al charge is gAi = 1-8, thereby producing charge neutrality. For the 
cut-off radius tai, the value tai = 6A is chosen. The parameter d^i is set to 2A^ for qAi{r) 
in the Al-0 interactions and to 1.47 for QAiif) in the Al-Al and the Al-Si interactions. 
Finally, we have chosen Cai = 0.0653609 for the Al-0 interactions and Cai = 0.0637977 for 
the Al-Al and Al-Si interactions. The parameters for c/ai and Cai are slightly different from 
those reported in Winkler et al. (2004). However, these slight differences are not relevant 
with respect to the resulting structural and dynamic properties. 

From Eq. it becomes obvious that at small distances the potential between the Al 
(or Si) and the O atoms goes to minus infinity (since the coefficients Ca/s are positive), 
i.e. it becomes unphysical. Therefore we have modified the potential at short distances by 
substituting it by a polynomial continuation that makes the potential repulsive at very short 
distances: 

(pao{r) = ai,a + a2,a r + ^a3,„(r - a^^o,)^ r < rf^^ (5) 

with a = Si, Al, O. The values for oi q,, 02,0, 03,0, 04,0, and the cut-off radius rP"'^ are 
listed in Table HI One may wonder why the poljTiomial terms in Eq. (0) are appropriate to 
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describe the interactions at small distances accurately. But even at the highest considered 
temperature, T = 6000 K, only a negligible number of ion pairs approach each other at 
distances r < rP"'^. Thus, the potentials as given by Eq. are only of technical importance 
and do not affect at all the physical properties in the temperature range considered in this 
study. 

Having described in detail the model potential, we give now the main details of 
the simulation runs. Molecular dynamics simulations were done for the two composi- 
tions 29(Al203)197(Si02) (with about 13mol% AI2O3) and 65(Al203)73(Si02) (with about 
47mol% AI2O3). In the following, we will abbreviate these compositions as A29S197 and 
A65S73, respectively. Both systems were simulated at two different mass densities p to 
study how pressure affects structural and dynamic properties. For A29S197, the densities 
p = 2.29 g/cm^ and p = 2.4g/cm^ were chosen, which are slightly below and above the 
experimental value at room temperature, p = 2.36 g/cm^ (Morikawa et ai, 1982). In the 
case of A65S73, the densities p = 2.35 g/cm'^ and p = 2.65 g/cm^ were chosen, both of which 
are below the experimental value p = 2.74g/cm^ (Morikawa et ai, 1982). In the following, 
we refer to runs at low and high density by the abbreviations LD and HD, respectively. 

The simulated systems consist of 2208 and 2176 atoms for A29S197 and A65S73, re- 
spectively, that sit in a cubic simulation box with periodic boundary conditions in all three 
Cartesian directions. The equations of motion were integrated with the velocity form of the 
Verlet algorithm, and the Coulombic contributions to the potential and the forces were calcu- 
lated via Ewald summation (Binder et ai, 2004). The time step of the integration was 1.6 fs. 
For each of the four systems (i.e. A29S197 and A65S73 at two different densities each), 13 
temperatures were considered in the interval 6000 K > T > 2000 K (the lowest temperatures 
were 2470 K and 2390 K for the LD and HD runs of A29S197, respectively, and 2190 K and 
2060 K for the LD and HD runs of A65S73, respectively). The temperature of the systems 
was controlled by coupling them to a stochastic heat bath, i.e. by periodically substituting 
the velocities of the particles with the ones from a Maxwell-Boltzmann distribution with 
the desired temperature. This thermostat has been first proposed by Andersen (1980) who 
has shown also that it generates a canonical distribution in phase space. After the system 
was equilibrated at the target temperature, we continued the run in the microcanonical 
ensemble, i.e. the heat bath was switched off. In order to improve the statistics we have 
done eight independent runs at each temperature. At the lowest temperatures, equilibration 
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runs between 23.4 ns and 42.3 ns real time were done, followed by production runs of the 
same length. The longest simulations were done for the HD A29S197 system at T = 2390 K, 
where the total simulation time was 2 x 42.3 ns (equilibration+production) x 8 = 676.8 ns, 
corresponding to about 416 million time steps. In addition to the simulations of fully equi- 
librated melts, we produced glass structures by quenching the systems from the melt to 
300 K, followed by relaxation runs of 10,000 time steps at that temperature. Note that the 
total simulation time of the present study was 48.2 CPU years on an IBM Regatta at the 
NIC Jiilich. 

Fig. HI displays the temperature dependence of the pressure p for the different systems. In 
the case of the HD samples, p{T) exhibits a minimum which moves to lower temperatures 
with increasing AI2O3 concentration, namely from about 4300 K in A29S197 to about 2400 K 
in A65S73. This feature is less pronounced in the LD samples, and it is even absent in the 
LD sample of A65S73. In agreement with density measurements of aluminium silicate melts 
around 2000 K, the pressure of the A29S197 model exhibits a relatively weak dependence 
on temperature. From Fig. one can infer that around 2000 K, the density at ambient 
pressure is about 2.3g/cm^. This value is in good agreement with the experimental value 
of 2.32 g/cm^ at 2000 K for an aluminium silicate melt with 14.82 mol% AI2O3 (Aksay et 
ai, 1979). However, for the Al203-Si02 melt with 47mol% AI2O3 (which is similar to our 
A65S73 model), Aksay et al. (1979) obtained a value of 2.62 g/cm^ at 2200 K. At this tem- 
perature, one can estimate an ambient pressure density around 2.45 g/cm^ for our simulation 
model and thus, in the case of the A65S73 model, the density is about 10% smaller than in 
the corresponding real system. 

III. RESULTS 

A. Structural properties 

In this section, structural properties of the models of A29S197 and A65S73 are investi- 
gated. The central issue is to elucidate the interplay between the local chemical ordering 
around the Al atoms and the IRO that we have recently encountered in a simulation study 
of AS2 (Winkler et a/., 2004). 

As we have already mentioned in the Introduction, the local structure of amorphous alu- 
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minium silicates has been investigated by various experimental techniques. In particular, 
NMR experiments yield detailed information about the Al-0 coordination and the occur- 
rence of triclusters. All these experiments have been done at room temperature, i.e. well 
below the experimental glass transition temperature Tg. Therefore, a direct comparison 
of, e.g., NMR results with those from a computer simulation is not very meaningful, since 
the glass structures from a MD simulation have a very different thermal history than the 
experimental ones (which is due to the much shorter time scale that is accessible in a simu- 
lation). All one can do is to extrapolate the properties of fully-equilibrated samples at high 
temperatures (in this study above 2000 K) to temperatures at which the system would fall 
out of equilibrium on the typical experimental time scale. 

Keeping this issue in mind, we discuss now the temperature dependence of the Al-0 
coordination in the A29S197 and A65S73 models for different coordination numbers zai-o- 
The coordination number -Zai-o is defined as the number of O atoms surrounding an alu- 
minium atom within a distance r < rcut- For rcut, we have chosen 2.32 A, which corresponds 
to the location of the first minimum in the partial pair correlation function of the Al-0 
correlations. Note that rcut is significantly larger than the mean distance between an Al and 
an O atom which is tai-o = 1-67 A for our model system. But this difference between r^ut 
amd tai-o is normal for a melt structure where the relatively large width of the first peaks 
in the pair correlation function is due to structural disorder and the thermal motion of the 
atoms. 

Fig. El shows the coordination number distribution Pai-o(^) in a semi-logarithmic plot 
as a function of inverse temperature for zm-o = 3, 4, 5. We have not included the results for 
zj^i_o = 6 since six-fold coordinated aluminium atoms are quite rare, i.e. not more than 2% 
of them are found in each of the considered systems. From our data it is hard to extrapolate 
to lower temperatures since none of the curves can be described by an Arrhenius law (a 
straight line in Fig. |2j over a large temperature range. It might even be that, e.g., the curves 
for z = 4 exhibit a maximum at large values of 1/T, or those for z = 5 a minimum. At least 
it seems that three-fold coordinated Al atoms disappear at low temperature whereas one 
may expect a significant fraction of five-fold coordinated Al atoms also at low temperature, 
in particular in the HD A65S73 system. However, Fig. |21 shows that for A29S197, as well as 
for A65S73, most of the Al atoms are four-fold coordinated by O atoms. 

Since five- and six-fold coordinated Al atoms are relatively rare, we expect a large num- 
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ber of triclusters of the form 0-3(Si,Al) in order to yield charge neutrahty in the local 
environment of Al atoms. And indeed, as can be seen in Fig. EK? the number of triclusters 
0-3(Si,Al) per Al atom is of the order of one in the whole temperature range that is con- 
sidered for the different systems. The minimum in the curves for A29S197 around 3000 K 
can be easily explained by a closer inspection of the data. At high temperature, triclusters 
with only one or no Al atom are the most frequent ones. But their number strongly de- 
creases with decreasing temperature, whereas the number of triclusters with two or three Al 
atoms increases with decreasing temperature, thus leading to the minima in Fig. E^. Fig. 
displays the number of two-fold rings (i.e. edge-sharing tetrahedra) as a function of temper- 
ature. Clearly, there are about 0.4 two-fold rings per Al atom at low temperature. Most of 
the two-fold rings contain two Al atoms (about 60-80%) or one Al and one Si atom (about 
20-40%). The fraction of edge-sharing Si04 tetrahedra is for all the considered systems of 
the order of 1% or smaller (of course, this holds only at low temperatures). The O atoms 
of the two-fold rings are most likely also triclusters, in agreement with our previous study 
(Winkler et a/., 2004 and Tossell and Horbach, 2005). 

In our recent work on AS2, we have demonstrated that the different chemical ordering 
of aluminium and silicon leads to a microphase separation on intermediate length scales of 
the order of 10-12 A whereby an Al-0 rich network percolates through the Si-0 network. 
These structural correlations are reflected by a prepeak in partial static structure factors 
at a wavenumber around q = 0.5 A~^. Appropriate quantities to study the latter feature 
are the static concentration-concentration structure factors Sc^cail)- In the following, we 
define these quantities and we discuss their behavior in the case of the aluminium silicates 
considered here. 

For the definition of the Sc^c^iQ) we follow a paper by Bletry (1976) that generalizes the 
Bhatia-Thornton formalism for two-component systems to arbitrarily many components. 
Consider a mixture of n chemical species that contains a total number of = J2a=i 
particles, where Na denotes the number of particles of species a. The local number density 
in reciprocal space for particles of type a is defined as follows (Hansen, 1986): 

Pa(q) = IIexp(iq- Ffc) (6) 

k=l 

with q the wavevector and the position of the /c'th particle of type a. The partial static 
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structure factors are then given by (Hansen, 1986) 

Here we assume that the system is isotropic and thus, each of the Sapiq) depends only 
on the magnitude of the wavevector q. The functions Saisiq) are obviously symmetric, 
i.e. Sa/siq) = SjSaiq), and therefore, there are n(n+l)/2 independent partial structure factors. 
With this definition, Eq. (jZj), the Sajsiq) approach in the limit q ^ oo the mean concentration 
Xa = Na/N for a = (3 and zero for a ^ (3. We show below that the concentration- 
concentration structure factors can be written as linear combinations of the partial structure 
factors Sap{q)- 

The local concentration variables for particles of type a are given by 

n 

Ca(q) = Pa(q) -XaY. pM ■ (8) 

/3=1 

The densities CQ,(q) express the local deviation from a homogeneous density distribution of 
particles of type a and thus, if pa is equal to Na/N for a = 1, n, the variables Cq vanish. 
The partial structure factors that correspond to the concentration densities c„(q) are defined 
in a similar way as the partial structure factors for the number densities: 

Sc^c,{q) = ^ (cMcfsi-q)) . (9) 

In this case Sc^c^iiq) = Sc^^cdq) also holds. Moreover, the functions Sc^c^iiq) obey the sum 
rule 

E^c.c,(g) = o (10) 

/3=1 

which follows directly from the definition, Eq. 0. 

In the case of n = 2 one has ci = — C2, and thus Sciaiq) = Sc2C2{q)- Furthermore, Eq. (fTnji 
yields Sec = Sc^^aiq) = —SciCc2{q)- This means that for n = 2 there is only one relevant 
structure factor Scc{q) for the concentration density correlations, and this quantity can be 
written as a linear combination of the partial structure factors as given by Eq. ((Tj), 

Scc{q) = xjSuiq) + xjS22{q) - 2xia;2S'i2(g) . (11) 
For g ^ oo this function approaches X1X2, which is expected for an ideal mixture. 
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Slightly more complicated is the case n = 3. Now the functions Sc^Cfsig) with a = P are 
related to the partial number density structure factors as follows: 

Sciciiq) = {x2 + xsf Su{q) -2xi{x2 + X3) [Suiq) + Si3{q)] 

+x?[^22(g) + 2523(g) + 533(g)] , (12) 

Sc2C2{q) = (a;i + X3)^5'22(g) - 2x2(0:1 + 0:3) [5*12 (g) + 5*23 (g)] 

+x2[5n(g) + 25i3(g) + 533(g)] , (13) 

Scscsiq) = {Xl + X2fS33{q) -2X3{XI+ X2)[Si3{q) + S23iq)] 

+x^[5n(g) + 25i2(g) + 522(g)] • (14) 
Each of the Sc^cp{q) with a ^ (3 can be written as a linear combination of the three 5c„cc,(g), 

\[Sc^cM - Sc^cAq) ~ Sc2C2{.q)] , (i5) 
\[Sc2C2{q) - Sc,cM - Scscaiq)] , (16) 
^[Sc^cM - Sc2C2{q) - Sccsiq)] ■ (i7) 

Thus, the latter functions do not contain any additional information, and so we consider 
only the correlation functions as given by Eqs. p2j) - (fT^ in the following. 

The three functions 5c^cQ(g) are shown in Fig. 0] for A29S197 and in Fig. El for A65S73. 
Their behavior is qualitatively very similar in the two systems: In both cases, the functions 
exhibit a pronounced peak around gi = 2.72 A~^. This peak indicates the chemical ordering 
on local length scales: Between nearest cation neighbors there is always a shell of oxygen 
atoms, i.e. strong fluctuations occur on length scales of the order of the nearest cation- 
oxygen distance (note that 27r/gi = 2.3 A corresponds approximately to the location of 
the first minima in the partial pair correlation function of Si-0 and Al-0 correlations). 
Apart from the peak at gi, in Sc^^cAiiq) ^csicsiiq), a pronounced prepeak is found around 
gp = 0.5 A^\ whereas in ScQcoiq), such a prepeak is absent. This behavior of the different 
concentration-concentration functions is similar to our recent result for AS2. There, we 
have related the prepeak to a microsegregation into Al-0 rich and Si-0 rich regions on the 
length scale of about 1 nm. Indeed, in the snapshots of Figs. E] and El one can clearly identify 
Al-0 rich percolating regions in the Si-0 network. It is remarkable that the location of the 
prepeak at gp seems to depend neither on the Al-0 concentration nor on the pressure (at 
least for the pressures that are observed at the two densities considered for each system). 
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Sc,c2{q) = 

Scicaiq) = 

Sc2cz{q) = 



This finding is contrary to the behavior of structural features on local length scales where 
we have found a strong dependence on temperature, pressure, and AI2O3 concentration (see 
above). 

The determination of Sc^cai^) from scattering experiments is very difficult since it requires 
the measurement of partial structure factors Sapiq). For silicates, this has been possible only 
for a few examples, e.g. a calcium silicate glass (Gaskell et ai, 1992). In the case of aluminium 
silicate glasses, total structure factors have been measured for different compositions using 
X-ray scattering by Morikawa et al. (1982) and by Okuno et al. (2005). Their results can 
be directly compared to simulation data. To this end, an X-ray scattering structure factor 
Sy^{q) has to be calculated from the partial structure factors, Sa/3{q), by weighting them 
with X-ray form factors, 

^x(g) = ^ E fa{s)Ms)S^f,{s) (18) 

with a, f3 = {Si, Al, O}. The form factors /a(s) depend on the wavenumber q via s = q/ (47r). 
Note that the fa{s) are taken from the literature (International Tables, 1974). Fig. IHlshows 
the "reduced" X-ray structure factor q{Sx{q) — 1) in comparison to experimental results. 
As we recognize from the figures a good agreement is obtained between simulation and 
experiment both for A29S197 and A65S73. For the latter system, two experimental data 
sets are available and, as we see in Fig. Eb, the simulation curve is closer to the more recent 
result of Okuno et al. (2005). 

As we can infer from Fig. no prepeak is visible around 0.5 A~^ in the reduced X-ray 
structure factor. Note that this is also the case for Sx{q) itself. This is due to the fact 
that the main contribution to Sx{q) comes from S'o-o(Q') which does not exhibit a prepeak 
at 0.5 A~^. Also in neutron scattering experiments, one has no access to partial structure 
factors for aluminium silicates (due to the lack of appropriate isotopes). Thus, the accessible 
total structure factor is dominated by the 0-0 correlations and will hardly show a prepeak. 
It remains a challenge to the experimentalists to verify the presence of the latter prepeak in 
aluminium silicates. 
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B. Dynamical properties 



We have seen that the considered aluminium sihcates exhibit a chemical ordering on 
length scales of about 1 nm which can be described by a microphase separation into an Al 
rich and a Si rich network structure. In our recent study of AS2 (Winkler et al., 2004), we 
have shown that the addition of AI2O3 to a silica melt enables a much faster self diffusion of 
all components compared to pure Si02. This is also the case for A29S197 and A65S73 which 
we consider here. In Fig. [7^, the selfdiffusion constants of our systems (for the LD samples) 
are displayed on a semilogarithmic scale as a function of inverse temperature. Clearly, the 
diffusion constants increase with increasing AI2O3 concentration (note that the diffusion 
constants of AS2 fall between the curves shown in Fig. [7^, as expected; see Winkler et a/., 
2004). In Fig. [7|d, the oxygen diffusion constants Dq are shown, now also with the results 
for the HD systems. We can infer from this figure that an anomalous behavior of Dq with 
respect to pressure emerges: The diffusion becomes faster with increasing pressure. This 
anomaly is well-known in many different "simple" network forming glasses such as Si02, 
H2O, etc. (see, e.g.. Shell et a/., 2002 and references therein). For Si02, this anomaly has been 
related to structural changes that occur with increasing pressure (Angell et al., 1982; Kubicki 
and Lasaga, 1988; Rustad et al., 1990; Barrat et al., 1997): The network structure changes 
gradually from a four-fold coordination of silicon atoms at low pressure to an imperfect 
five-fold coordination at intermediate pressure and eventually to a six-fold coordination at 
very high pressure (note that in the latter regime the selfdiffusion constants decrease with 
increasing pressure). Intuitively, the anomalous diffusion can be understood as follows: In 
the tetrahedral network structure there are not as many possibilities for the formation of 
local defects as in a higher-coordinated network. Since these defects facilitate diffusion in a 
network structure (see, e.g., Horbach and Kob, 1999), the particles in a higher-coordinated 
network tend to exhibit a faster diffusion than in a tetrahedral network. 

But what are the defects that facilitate the diffusion in the Al203-Si02 networks? Im- 
portant structural units might be triclusters and five-fold coordinated sihcon atoms, both 
of which are more frequent in the HD systems (see above). Thus, they might be responsible 
for the higher diffusion constants, as compared to the LD systems. 

As we see in Fig. [7^, silicon is the slowest component and aluminium the fastest com- 
ponent of the diffusing species. In order to quantify the temperature dependence of the 
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different diffusion constants, we show in Fig. |Slthe ratios Dsi/Do and Dm/Dq as a function 
of temperature for tfie different systems. It is remarkable that Da,i/ Dq depends only weakly 
on density and composition (note that also D^i/Do for BKS-Si02 would fall roughly on 
top of the corresponding curves in Fig. |H1 see Winkler et ai, 2004). This shows that the 
diffusion of silicon and oxygen is intimately connected with each other. On the other hand, 
Dai/ Do nearly approaches the constant value of one in A65S73 in the considered temper- 
ature range (at least for the HD system). This might be due to the fact that in the AI2O3 
rich system, A65S73, oxygen diffusion steps occur mainly in the vicinity of Al rich regions, 
in particular near triclusters, whereas in A29S197, also the O diffusion in the slow Si rich 
regions is important. This issue has to be clarified in future studies. 



IV. SUMMARY 



Large scale molecular dynamics computer simulations were used to study the structure 
and diffusion dynamics of the aluminium silicates A29S197 and A65S73. The microscopic 
interactions between the ions were described by a simple pair potential proposed by Kramer 
et al. (1991). In this work, we have demonstrated that this model potential yields good 
agreement with available experimental data such as the mass density and the X-ray structure 
factor. However, recent NMR experiments (Sen and Youngman, 2004; Stebbins et al, 2005) 
observed a significant number of six-fold coordinated Al atoms already at small AI2O3 
concentrations. By contrast, our simulation model predicts a vanishing number of AlOg 
units in the glass structure. This might be realistic for A29S197 but not for the Al-rich 
A65S73 system. In the case of five-fold coordinated Al atoms, an extrapolation from the 
melt data to the experimental glass transition temperature is not easy (see above). However, 
as a typical local structure, we find tricluster O atoms in conjunction with two-fold rings. 
Recent Hartree-Fock calculations that used small clusters from our recent MD simulation 
of AS2, found evidence that the combination of triclusters with two-fold rings is indeed 
realistic (Tossell and Horbach, 2005). It should be emphasized that this finding refers only 
to binary alumino silicates. The local structure may be very different in ternary alumino 
silicates that contain also charge-balancing cations such as calcium or alkali ion species. 

The different local ordering of Al and Si atoms leads to a microsegregation on length 
scales of about 1 nm. This structural feature is reflected in a prepeak in static concentration- 
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concentration correlation functions (and also in partial static structure factors). The location 
of the prepeak is around 0.5 A~^, independent of AI2O3 concentration, pressure, and the 
details of the local structure in the considered systems. This remarkable result shows that 
according to our simulation model, the formation of Al rich channels that percolate through 
the Si-0 network is a central feature of the chemical ordering in the system Al203-Si02. This 
is a prediction of our simulation that could be tested by scattering experiments. Prepeaks in 
other silicates, that have a similar origin as the one found in our aluminium silicate models, 
have been successfully identified in recent neutron scattering experiments (see Gaskell et al, 
1991; Meyer et al, 2002; Meyer et al, 2004; Kargl and Meyer, 2005). 

Our analysis of the diffusion dynamics suggests that triclusters are important for the 
mass transport in aluminium silicates. Moreover, wc have found subtle differences in the 
dependence of the three diffusion constants on temperature (e.g. the ratio Dai/Dq is essen- 
tially one for the HD A65S73 system, while it increases in the case of the A29S197 systems). 
It would be interesting to see whether these features can be understood in the framework of 
mode coupling theory. Work in this direction is in progress. 
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VI. LIST OF CAPTIONS 



Fig. [TJ Pressure as a function of temperature at the indicated compositions and densi- 
ties. Note that the points at 300 K are far below the glass transition temperature of the 
simulations, which is around 2000 K (see text). 

Fig. 121 Temperature dependence of Pai-o(^) for z = 3, 4, 5 as indicated, for a) A29S197 
and b) A65S73. 

Fig. 01 This plot shows the amount of triclusters, 0-3(Si,Al), and 2-fold rings as a 
function of temperature: a) Number of triclusters and b) number of 2-fold rings per Al 
atom for the different systems as indicated. 

Fig. m Structure factors Sc^cai^) fo^' A29S197 at the temperature T = 2470 K and the 
different densities as indicated, a) S'cgicsi(q'), b) Scqco{(i), and c) Sc^jc^i (g) [see Eqs. 
for the definition of Sc^caiQ)]- Also shown is a snapshot at T = 300 K which illustrates 
the intermediate range order as reflected in the prepeak around q = 0.5 A~^. The large 
white spheres are the silicon atoms, the large black spheres are the aluminium atoms, and 
the small black spheres are the oxygen atoms. Note that the size of the spheres does not 
correspond to the actual size of the atoms. 

Fig. El The same as in Fig. |3]but now for A65S73. The temperatures are T = 2190 K 
and T = 2060 K for the low and the high density, respectively. 

Fig. ini Reduced X-ray scattering structure factor q{Sx{(l) — 1) as calculated from the 
simulation using Eq. (fTH|) in comparison to experimental results by Morikawa et al. (1982) 
and Okuno et al. (2005) (dashed lines), for a) A29S197 and b) A65S73. 

Fig. [71 Arrhenius plots of the selfdiffusion constants D^, a) for Si, Al, and O for the 
different compositions at low densities and b) for O for all considered systems as indicated. 

Fig. ISl Temperature dependence of the ratios of the diffusion constants. Da/ Dp (a,/? = 
Si, Al, O), for the different systems as indicated. 
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VII. LIST OF TABLES 



a-0 


ai,a [eV] 


a2,a [eV/A] 


a3,« [eV/A3] 


a4,a [A] 




Si-0 


-23.96027 


-2.85441 


-50.0 


1.41590 


1.276 


Al-0 


-87.62405 


26.72474 


-3.0 


4.49012 


1.28 





-30590.38 


90.38499 


-0.1 


97.25877 


1.9547 



TABLE I: Parameters for the polynomial continuation at small r that prevents the potentials from 
going to minus infinity for r ^ 0. 
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